#############################################################################
# Modeling Issue Competence Over Time
# Last Updated: 2025-12
# R version 4.5.1 (2025-06-13) -- "Great Square Root"
# Input:
#         1. full_base_rev_v2.Rdata (Model 1 in Table 1)
#         2. full_base_pid_rev_v2.Rdata (Model 2 in Table 1)
#         3. full_topic_rev_v2.Rdata (Model 3 in Table 1)
#         4. full_topic_pid_rev_v2.Rdata (Model 4 in Table 1)
#         5. full_topic_biannual_interact_pid_rev_v2.Rdata (Dynamic Issue Competence model; Table A2)
#         6. issue_static_list.rds: Model prediction by Model 4 (will be created during the code)
#         7. issue_trend_list.rds: Model predictions by dynamic issue competence model (will be created during the code)
############################################################################

rm(list = ls())

con <- file("replication_log.txt")
sink(con, append=TRUE)
sink(con, append=TRUE, type="message")

#-----------------------------
# Load packages and data
#-----------------------------
library(brms)
library(xtable)
library(stringr)
library(ggplot2)
library(tidyr)
library(dplyr)
library(tibble)
library(dotwhisker)
library(tidytext)

########################################################################
#Static Issue Competence
########################################################################
load("full_base_rev_v2.Rdata") #Model 1
load("full_base_pid_rev_v2.Rdata") #Model 2
load("full_topic_rev_v2.Rdata") #Model 3
load("full_topic_pid_rev_v2.Rdata") #Model 4

#Report Results (Table 1)
  #Model 1
  dat <- full_base$data
  model <- full_base
  round(summary(model)$fixed, 3)
  my_results<-round(summary(model)$fixed, 3)
  write.table(my_results, file = "table1_model1.txt", sep = "\t", quote = FALSE, col.names = NA)
  
  #Model 2
  dat1 <- full_base_pid$data
  model <- full_base_pid
  round(summary(model)$fixed, 3)
  my_results<-round(summary(model)$fixed, 3)
  write.table(my_results, file = "table1_model2.txt", sep = "\t", quote = FALSE, col.names = NA)
  
  #Model 3
  dat <- full_topic$data
  model <- full_topic
  round(summary(model)$fixed, 3)
  round(summary(model)$random$cap_full, 3)
  my_results1<-round(summary(model)$fixed, 3) 
  my_results2<-round(summary(model)$random$cap_full, 3) 
  my_results<-rbind(my_results1, my_results2)
  write.table(my_results, file = "table1_model3.txt", sep = "\t", quote = FALSE, col.names = NA)
  
    #Random effect (Table A3; Model 3)
    re <- ranef(model, summary = TRUE, probs = c(0.025, 0.975))$cap_full
    round(re, 3)
    my_results<-round(re, 3)
    write.table(my_results, file = "tableA3_1.txt", sep = "\t", quote = FALSE, col.names = NA)
  
  #Model 4
  dat <- full_topic_pid$data
  model <- full_topic_pid
  round(summary(model)$fixed, 3)
  round(summary(model)$random$cap_full, 3)
  my_results1<-round(summary(model)$fixed, 3) # Estimate and Error
  my_results2<-round(summary(model)$random$cap_full, 3) #Group-level Issue Estimate and Error
  my_results<-rbind(my_results1, my_results2)
  write.table(my_results, file = "table1_model4.txt", sep = "\t", quote = FALSE, col.names = NA)
  
  #Random effect (Table A3; Model 4)
  re <- ranef(model, summary = TRUE, probs = c(0.025, 0.975))$cap_full
  round(re, 3)
  my_results<-round(re, 3)
  write.table(my_results, file = "tableA3_2.txt", sep = "\t", quote = FALSE, col.names = NA)
  
  #Visualization of Model 4 (Figure A1)
  cols <- c(
    "CDU/CSU"   = rgb(0, 0, 0),         # CDU/CSU
    "SPD"       = rgb(227/255, 0, 15/255),
    "Grüne"     = rgb(70/255, 150/255, 43/255),
    "FDP"       = rgb(255/255, 222/255, 0),
    "Die Linke" = rgb(207/255, 52/255, 118/255),
    "AfD"       = rgb(0, 156/255, 222/255),
    "None"      = "grey60"
  )
  
  fx <- as.data.frame(fixef(model, probs = c(.025, .975))) %>%
    rownames_to_column("par")
  
  pid_fx <- fx %>%
    filter(str_detect(par, "_pid")) %>%
    separate(par, into = c("resp_code", "pid_code"), sep = "_", remove = FALSE)
  
  party_levels <- c("SPD","Grüne","FDP","Die Linke","AfD")
  pid_levels   <- c("CDU/CSU","SPD","Grüne","FDP","Die Linke","AfD","None")
  
  pid_fx <- pid_fx %>%
    mutate(
      party = recode(resp_code,
                     mu02spd="SPD", mu03gruene="Grüne", mu04fdp="FDP",
                     mu05linke="Die Linke", mu06afd="AfD"),
      pid_label = recode(pid_code,
                         pid01cdu="CDU/CSU", pid02spd="SPD", pid03gruene="Grüne",
                         pid04fdp="FDP", pid05linke="Die Linke",
                         pid06afd="AfD", pid07none="None"),
    
      OR = exp(Estimate), lo = exp(Q2.5), hi = exp(Q97.5),
      party    = factor(party, levels = party_levels),
      pid_label= factor(pid_label, levels = pid_levels)
    )
  
  p1 <- ggplot(pid_fx, aes(x = OR, y = pid_label, color = pid_label)) +
    geom_vline(xintercept = 1, linetype = "dashed") +
    geom_point(size = 2) +
    geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.2) +
    facet_wrap(~ party, ncol = 3, scales = "free_y") +
    scale_x_log10() +
    scale_color_manual(values = cols) +
    labs(
      x = "Odds ratio", y = "Respondent PID",
      color = "Respondent PID"
    ) +
    theme_minimal() +
    theme(legend.position = "none")
  p1
  
  ggsave(
    filename = "group_fixed_effect.jpeg",
    plot = p1,
    width = 8,
    height = 6 
  )
  
# Model Predictions: Model predictions: Static issue competences along with 95 percent credible interval
  # Figure 1 & Table A4
  
  pids_table  <- table(dat$pid)
  pid_list    <- c("01cdu","02spd","03gruene","04fdp","05linke","06afd","07none")
  party_names <- c("CDU/CSU","SPD","Grüne","FDP","Die Linke","AfD","None")
  party_names6 <- c("CDU/CSU","SPD","Grüne","FDP","Die Linke","AfD")
  
  w_all <- as.numeric(pids_table[pid_list])
  w_all <- w_all / sum(w_all)
  
#   out_list <- list()
#   for (issue in unique(dat$cap_full)) {
#     
#     pred_dat <- expand.grid(
#       cap_full = issue,
#       pid      = pid_list
#     )
#     
#     # posterior epred: draws × G × K
#     pred <- posterior_epred(model, newdata = pred_dat, allow_new_levels = TRUE)
#     D <- dim(pred)[1] #4000 predictions
#     G <- dim(pred)[2] 
#     K <- dim(pred)[3]
#     
#     mixed <- matrix(NA_real_, nrow = D, ncol = K)
#     for (d in seq_len(D)) {
#       mixed[d, ] <- drop(crossprod(w_all, pred[d, , ]))
#     }
#     
#     out_mean <- colMeans(mixed)
#     out_lo   <- apply(mixed, 2, quantile, 0.025)
#     out_hi   <- apply(mixed, 2, quantile, 0.975)
#     
#     out_mean <- matrix(out_mean, nrow = K, ncol = 1, dimnames = list(party_names6, NULL))
#     out_lo   <- matrix(out_lo,   nrow = K, ncol = 1, dimnames = list(party_names6, NULL))
#     out_hi   <- matrix(out_hi,   nrow = K, ncol = 1, dimnames = list(party_names6, NULL))
#     
#     out_list[[issue]] <- list(
#         mean = out_mean,
#         lo   = out_lo,
#         hi   = out_hi
#       )
#   }
#   
#   saveRDS(out_list, "issue_static_list.rds")
  
  out_list <- readRDS("issue_static_list.rds")
  unique(dat$cap_full)
  
  #Table A4
  round(out_list[["01_macroeconomics"]]$mean,3)
  round(out_list[["01_macroeconomics"]]$lo,3)
  round(out_list[["01_macroeconomics"]]$hi,3)
  
    sorted_issues <- sort(unique(dat$cap_full))
    final_table <- NULL
    for (issue_name in sorted_issues) {
      
      if (!issue_name %in% names(out_list)) next 
      
      item <- out_list[[issue_name]]
      val_mean <- item$mean
      val_lo   <- item$lo
      val_hi   <- item$hi
      
      party_names <- rownames(val_mean)
      cell_content <- sprintf("%.3f\n[%.3f, %.3f]", val_mean, val_lo, val_hi)
      
      if (is.null(final_table)) {
        final_table <- data.frame(t(cell_content))
        colnames(final_table) <- party_names
        
      } else {
        new_row <- data.frame(t(cell_content))
        colnames(new_row) <- colnames(final_table)
        final_table <- rbind(final_table, new_row)
      }
    }
    
    clean_names <- gsub("^[0-9]+_", "", sorted_issues)
    rownames(final_table) <- clean_names
    
    write.table(final_table, "Table_A4_result.txt", sep="\t", quote=FALSE, col.names=NA)
    
  #Visualize - Figure 1
  out_df <- purrr::imap_dfr(out_list, function(issue_list, issue_name) {
    # Extract vectors
    mean_df <- as.data.frame(issue_list$mean)
    lo_df   <- as.data.frame(issue_list$lo)
    hi_df   <- as.data.frame(issue_list$hi)
    
    # Get party names
    parties <- rownames(mean_df)
    
    # Combine into a tidy tibble
    tibble(
      issue = issue_name,
      party = parties,
      mean  = mean_df[[1]],
      lo    = lo_df[[1]],
      hi    = hi_df[[1]]
    )
  })
  
  party_order <- c("CDU/CSU","SPD","Grüne","FDP","Die Linke","AfD")
  out_df <-out_df %>% dplyr::mutate(party=factor(party, levels = party_order))
  
  out_df <- out_df %>%
    mutate(
      issue = case_when(
        issue == "01_macroeconomics" ~ "Macroeconomics",
        issue == "02_civilrights"    ~ "Civil Rights",
        issue == "03_health"         ~ "Health",
        issue == "05_labor"          ~ "Labor",
        issue == "06_education"      ~ "Education",
        issue == "07_environment"    ~ "Environment",
        issue == "08_energy"         ~ "Energy",
        issue == "09_immigration"    ~ "Immigration",
        issue == "10_transportation" ~ "Transportation",
        issue == "12_lawcrime"       ~ "Law and Crime",
        issue == "13_welfare"        ~ "Welfare",
        issue == "14_housing"        ~ "Housing",
        issue == "15_commerce"       ~ "Commerce",
        issue == "16_defense"        ~ "Defense",
        issue == "19_international"  ~ "International",
        issue == "20_government"     ~ "Government",
        TRUE                         ~ issue
      )
    )
  
  cols <- c(
    rgb(0, 0, 0),
    rgb(227/255, 0/255, 15/255),
    rgb(70/255, 150/255, 43/255),
    rgb(255/255, 222/255, 0/255),
    rgb(207/255, 52/255, 118/255),
    rgb(0/255, 156/255, 222/255)
  )
  names(cols) <- c("CDU/CSU","SPD","Grüne","FDP","Die Linke","AfD")
  
  p <- out_df %>%
    dplyr::mutate(issue_in_party = tidytext::reorder_within(issue, mean, party)) %>%
    ggplot(aes(x = mean, y = issue_in_party, color = party)) +
    geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0.2, size = 0.5) +  
    geom_point(size = 1.6) +                                              
    facet_wrap(~ party, scales = "free_y", ncol = 2) +
    tidytext::scale_y_reordered() +
    scale_color_manual(values = cols) +                            
    labs(x = "Predicted Proportion", y = "Issue") +
    theme_minimal(base_size = 13)+
    theme(legend.position = "none")
  p
  ggsave(p, filename = "static_issue_competence.jpeg", width = 8, height = 10, dpi = 300)
  
########################################################################
#Dynamic Issue Competence
########################################################################

load("full_topic_biannual_interact_pid_rev_v2.Rdata") #Table A2
  
  #Report Results (Table A2)
  dat <- full_topic_biannual_interact_pid$data
  model <- full_topic_biannual_interact_pid
  summary(model)
  round(summary(model)$fixed, 3)
  round(summary(model)$random$cap_full, 3) #Topic-level
  round(summary(model)$random$halbjahr, 3) #Half-year
  round(summary(model)$random$"cap_full:halbjahr", 3) #Topic x Half-year
  my_results1<-round(summary(model)$fixed, 3)
  my_results2<-round(summary(model)$random$cap_full, 3) #Topic-level
  my_results3<-round(summary(model)$random$halbjahr, 3) #Half-year
  my_results4<-round(summary(model)$random$"cap_full:halbjahr", 3) #Topic x Half-year
  my_results<-rbind(my_results1, my_results2,my_results3,my_results4)
  write.table(my_results, file = "tableA2.txt", sep = "\t", quote = FALSE, col.names = NA)
  
  #Random effect
  re <- ranef(model, summary = TRUE, probs = c(0.025, 0.975))$cap_full
  round(re, 3)
  
  ## Summary of the raw data
    table(dat$cap_full)  # 16 issues
    
    # Party x Topic
    table(dat$cap_full[dat$pid == "01cdu"])
    table(dat$cap_full[dat$pid == "02spd"])
    table(dat$cap_full[dat$pid == "03gruene"])
    table(dat$cap_full[dat$pid == "04fdp"])
    table(dat$cap_full[dat$pid == "05linke"])
    table(dat$cap_full[dat$pid == "06afd"])
    table(dat$cap_full[dat$pid == "07none"])
    
    ## Year x PID
    xtable(table(dat$halbjahr, dat$pid))
  
  #-----------------------------
  # Issue-specific Predictions
  #-----------------------------
  pids_year <- table(dat$pid, dat$halbjahr)
  hj_labels <- colnames(pids_year)
  hj_years <- as.numeric(sub("_2", ".5", sub("_1", ".0", hj_labels)))
  
  party_names <- c("CDU/CSU", "SPD", "Grüne", "FDP", "Die Linke", "AfD")
  pid_list   <- c("01cdu", "02spd", "03gruene", "04fdp", "05linke", "06afd")
  pid_pred    <- c(pid_list, "07none")
  n_parties <- length(party_names)
  n_periods <- ncol(pids_year)
  
  # out_list <- list()
  # 
  # for (issue in unique(dat$cap_full)) {
  # 
  #   pred_dat <- expand.grid(
  #     cap_full = issue,
  #     halbjahr = hj_labels,
  #     pid = pid_pred
  #   )
  # 
  #   pred <- posterior_epred(model, newdata = pred_dat, allow_new_levels = TRUE)
  # 
  #   D <- dim(pred)[1] #Draw
  #   Nnew <- dim(pred)[2]
  #   K <- dim(pred)[3] #poutcome
  # 
  #   out_mean <- matrix(NA_real_, nrow = K, ncol = n_periods)
  #   out_lo   <- matrix(NA_real_, nrow = K, ncol = n_periods)
  #   out_hi   <- matrix(NA_real_, nrow = K, ncol = n_periods)
  # 
  #   for (i in seq_len(n_periods)) {
  #     this_hj <- hj_labels[i]
  #     rows_i  <- which(pred_dat$halbjahr == this_hj)
  # 
  #     w <- as.numeric(pids_year[pid_pred, i])
  #     w <- w / sum(w)
  # 
  #     mixed <- matrix(NA_real_, nrow = D, ncol = K)
  #     for (d in seq_len(D)) {
  #       mixed[d, ] <- drop(crossprod(w, pred[d, rows_i, ]))
  #     }
  # 
  #     out_mean[, i] <- colMeans(mixed)
  #     out_lo[, i]   <- apply(mixed, 2, quantile, 0.025)
  #     out_hi[, i]   <- apply(mixed, 2, quantile, 0.975)
  #   }
  # 
  #   rownames(out_mean) <- rownames(out_lo) <- rownames(out_hi) <- party_names
  #   colnames(out_mean) <- colnames(out_lo) <- colnames(out_hi) <- hj_labels
  # 
  #   out_list[[issue]] <- list(
  #     mean = out_mean,
  #     lo   = out_lo,
  #     hi   = out_hi
  #   )
  # }
  # 
  # saveRDS(out_list, "issue_trend_list.rds")

  # Visualization: Time panel (Figure 2)
  out_list <- readRDS("issue_trend_list.rds")
  
  alpha <- 0.1
  
  cols_alpha <- c(
    rgb(0, 0, 0, alpha),
    rgb(227/255, 0/255, 15/255, alpha),
    rgb(70/255, 150/255, 43/255, alpha),
    rgb(255/255, 222/255, 0/255, alpha),
    rgb(207/255, 52/255, 118/255, alpha),
    rgb(0/255, 156/255, 222/255, alpha)
  )
  names(cols_alpha) <- c("01cdu", "02spd", "03gruene", "04fdp", "05linke", "06afd")
  
  cols <- c(
    rgb(0, 0, 0),
    rgb(227/255, 0/255, 15/255),
    rgb(70/255, 150/255, 43/255),
    rgb(255/255, 222/255, 0/255),
    rgb(207/255, 52/255, 118/255),
    rgb(0/255, 156/255, 222/255)
  )
  names(cols) <- c("01cdu", "02spd", "03gruene", "04fdp", "05linke", "06afd")

  for (issue1 in names(out_list)) {
    
    if (issue1 %in% c("03_health", "09_immigration", "10_transportation")) {
      
      # Predicted values
      t_out_mean <- out_list[[issue1]]$mean
      t_out_lo   <- out_list[[issue1]]$lo
      t_out_hi   <- out_list[[issue1]]$hi
      
      rownames(t_out_mean) <- rownames(t_out_lo) <- rownames(t_out_hi) <- party_names
      
      # File name
      if (!dir.exists("issue_trends")) {
        dir.create("issue_trends")
      }
      jpeg(filename = paste0("issue_trends/", issue1, "_predicted_panel.jpeg"),
           width = 8, height = 6, units = "in", res = 300)
      
      par(mar = c(4, 4, 2, 1))  # margin: bottom, left, top, right
      
      # Empty plot
      plot(NA,
           xlim = range(hj_years),
           ylim = c(0, 1),
           xlab = "Year",
           ylab = "Predicted Proportion",
           #main = paste("Issue:", issue1),
           xaxt = "n", bty = "n")
      axis(1, at = seq(2009, 2023, by = 1))
      
      # Loop over parties
      for (j in seq_along(party_names)) {
        party <- party_names[j]
        pid <- pid_list[j]
        
        y_mean <- t_out_mean[party, ]
        y_lo   <- t_out_lo[party, ]
        y_hi   <- t_out_hi[party, ]
        
        if (party == "AfD") {
          y_mean[hj_years < 2013] <- NA
          y_lo[hj_years < 2013] <- NA
          y_hi[hj_years < 2013] <- NA
        }
        
        polygon(
          c(hj_years, rev(hj_years)),
          c(y_lo, rev(y_hi)),
          border = NA,
          col = cols_alpha[pid]
        )
        
        lines(hj_years, y_lo, col = 'darkgray', lwd = 1)  # lower bound
        lines(hj_years, y_hi, col = 'darkgray', lwd = 1)  # upper bound
        
        lines(hj_years, y_mean, col = cols[pid], lwd = 2)
        points(hj_years, y_mean, col = cols[pid], pch = 16)
      }
      
      legend(
        "topleft",
        legend = party_names,
        fill = cols[pid_list],     
        border = NA,               
        bty = "n"                 
      )
      
      dev.off()
    }
  }
  
  # Map issues to the single party (by pid code) you want to plot
  target_pairs <- c(
    "10_transportation" = "01cdu",   # CDU/CSU
    "14_housing"        = "02spd",   # SPD
    "13_welfare"        = "05linke", # Die Linke
    "09_immigration"    = "06afd"    # AfD
  )
  
  for (issue1 in names(out_list)) {
    
    # Only consider the four issues of interest
    if (issue1 %in% c("10_transportation", "14_housing", "13_welfare", "09_immigration")) {
      
      # Which party (index j) corresponds to this issue?
      pid_target <- unname(target_pairs[issue1])
      j <- match(pid_target, names(cols))    # index in your party vectors
      if (is.na(j)) next                     # safety guard if not found
      
      # Get predicted values
      t_out_mean <- out_list[[issue1]]$mean
      t_out_lo   <- out_list[[issue1]]$lo
      t_out_hi   <- out_list[[issue1]]$hi
      
      # Get raw values for this issue
      all_combinations <- expand.grid(
        pid = names(cols),      # your pid codes: "01cdu", ...
        halbjahr = hj_labels    # all half-year labels
      )
      
      dat_issue <- subset(dat, cap_full == issue1)
      
      dat_issue_counts <- dat_issue %>%
        count(pid, halbjahr, name = "count") %>%
        right_join(all_combinations, by = c("pid", "halbjahr")) %>%
        mutate(count = ifelse(is.na(count), 0, count))
      
      raw_table <- xtabs(count ~ pid + halbjahr, data = dat_issue_counts)
      raw_prop  <- prop.table(raw_table, margin = 2)
      
      raw_subset <- as.data.frame.matrix(raw_prop)
      raw_subset$pid <- rownames(raw_prop)
      rownames(raw_subset) <- NULL
      
      rownames(t_out_mean) <- rownames(t_out_lo) <- rownames(t_out_hi) <- party_names
      
      # Pull the target party info
      party_name <- party_names[j]
      pid_code   <- names(cols)[j]  # e.g., "01cdu"
      
      y_mean <- t_out_mean[party_name, ]
      y_lo   <- t_out_lo[party_name, ]
      y_hi   <- t_out_hi[party_name, ]
      
      if (pid_code == "06afd") {
        y_mean[hj_years < 2013] <- NA
        y_lo[hj_years < 2013]   <- NA
        y_hi[hj_years < 2013]   <- NA
      }
      
      # Clean filename and ensure folder
      party_filename <- gsub("/", "-", party_name)
      if (!dir.exists("comparison_graph")) dir.create("comparison_graph")
      
      jpeg(filename = paste0("comparison_graph/", issue1, "_", party_filename, "_compare.jpeg"),
           width = 8, height = 6, units = "in", res = 300)
      
      par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
      
      ## Left: Predicted panel
      plot(
        NA, xlim = range(hj_years), ylim = c(0, 1),
        ylab = "Predicted Proportion", xlab = "Year",
        bty = "n", xaxt = "n"
      )
      axis(1, at = seq(2009, 2023, by = 1))
      
      polygon(
        c(hj_years, rev(hj_years)),
        c(y_lo, rev(y_hi)),
        border = NA,
        col = adjustcolor(cols[j], alpha.f = 0.2)
      )
      points(hj_years, y_mean, col = cols[j], pch = 16)
      lines(hj_years, y_mean, col = cols[j])
      
      ## Right: Raw values panel
      plot(
        NA, xlim = range(hj_years), ylim = c(0, 1),
        ylab = "Raw Proportion", xlab = "Year",
        bty = "n", xaxt = "n"
      )
      axis(1, at = seq(2009, 2023, by = 1))
      
      if (pid_code %in% raw_subset$pid) {
        raw_row <- raw_subset %>% dplyr::filter(pid == pid_code) %>% dplyr::select(-pid)
        t_vector <- as.numeric(raw_row[1, ])
        if (pid_code == "06afd") t_vector[hj_years < 2013] <- NA
        
        points(hj_years, t_vector, col = cols[j], pch = 16)
        lines(hj_years, t_vector, col = cols[j])
      }
      
      dev.off()
    }
  }
  
  #Raw estimate
  raw_data <- data.frame()
  
  #For visualization
    cols_named <- c(
      "CDU/CSU"   = cols[["01cdu"]],
      "SPD"       = cols[["02spd"]],
      "Grüne"     = cols[["03gruene"]],
      "FDP"       = cols[["04fdp"]],
      "Die Linke" = cols[["05linke"]],
      "AfD"       = cols[["06afd"]]
    )
  
  #Combine raw data
  for (issue1 in names(out_list)) {
    
    if (issue1 != "Macro") {
     
      # Get raw values for this issue
      all_combinations <- expand.grid(
        pid = names(cols),  # your pid codes: "01cdu", ...
        halbjahr = hj_labels  # all half-year labels
      )
      all_combinations <- all_combinations %>% 
        mutate(hj_years = as.numeric(str_replace(str_replace(halbjahr, "_2$", ".5"), "_1$", ".0")))
      
      dat_issue <- subset(dat, cap_full == issue1)
      
      dat_issue_counts <- dat_issue %>%
        count(pid, halbjahr, name = "count") %>%
        right_join(all_combinations, by = c("pid", "halbjahr")) %>%
        mutate(count = ifelse(is.na(count), 0, count))
      
      total <- dat_issue_counts %>% group_by(halbjahr) %>% summarise(total=sum(count))
      dat_issue_counts <-dat_issue_counts %>% left_join(total, by = "halbjahr")
     
      dat_issue_counts <-dat_issue_counts %>% 
        mutate(prob = if_else(total==0, NA, count/total),
               issue = issue1,
               pid = case_when(
                 pid=="01cdu" ~ "CDU/CSU",
                 pid=="02spd" ~ "SPD",
                 pid=="03gruene" ~ "Grüne",
                 pid=="04fdp" ~ "FDP",
                 pid=="05linke" ~ "Die Linke",
                 pid=="06afd" ~ "AfD"
               ))
  
      raw_data<-rbind(raw_data, dat_issue_counts)
      }
  }
  #Afd
  raw_data <- raw_data %>% mutate(count = if_else(pid=="AfD" & hj_years<2013, NA_real_, count),
                                  prob = if_else(pid=="AfD" & hj_years<2013, NA_real_, prob))
  
  tab <- table(raw_data$count)
  sum(tab[names(tab) <= 10]) / sum(tab)
  median(raw_data$count, na.rm = TRUE)
  
  #Mean response for Linke on welfare
  linke <- raw_data %>% filter(pid=="Die Linke") %>% filter(issue=="13_welfare")
  mean(linke$count, na.rm = TRUE)

  #Mean response for afd on immigration
  afd <- raw_data %>% filter(pid=="AfD") %>% filter(issue=="09_immigration")
  mean(afd$count, na.rm = TRUE)
 
  #Visualize the number of raw data for each point (Figure C1)
  png(filename = "figure_c2.png", width = 800, height = 600)
  bp<-barplot(tab / sum(tab), xlab = "Count", ylab = "Proportion", xaxt="n")
  cats_num <- as.numeric(names(tab))
  ticks <- c(10, 100, 500, 1000)
  ticks_in_data <- ticks[ticks %in% cats_num]
  idx <- match(ticks_in_data, cats_num)
  axis(1, at = bp[idx], labels = ticks_in_data)
  dev.off()
  
  #Visualization: Comparison within issue (Figure C2)
    result <- raw_data %>%
      mutate(prob = ifelse(pid=="AfD" & hj_years < 2013, NA_real_, prob)) %>%
      filter(!is.na(prob)) %>%
      group_by(issue, halbjahr, hj_years) %>%
      slice_max(order_by = prob, n = 1, with_ties = TRUE) %>%
      ungroup() %>%
      dplyr::select(issue, halbjahr, hj_years, pid, prob, total)
    result
  
    #Find issue owner
    result <- result %>%
      group_by(issue) %>%
      summarise(winner = unique(pid), .groups = "drop") %>%
      distinct()
    
    for (issue1 in unique(raw_data$issue)) {
      
      if (issue1 != "10_transportation") next
      
      sub <-raw_data %>% filter(issue==issue1)
      
      party_series <- lapply(party_names, function(party) {
        s <- sub %>% filter(pid == party)
        m <- setNames(s$prob, s$hj_years)
        y <- vapply(hj_years, function(x) if (!is.na(m[as.character(x)])) m[as.character(x)] else NA_real_, 0.0)
        y[!names(y) %in% names(m)] <- NA_real_
        as.numeric(y)
      })
      names(party_series) <- party_names
      
      if (!dir.exists("raw_issue_trends")) {
        dir.create("raw_issue_trends")
      }
      jpeg(filename = paste0("raw_issue_trends/", issue1, "_raw_panel.jpeg"),
           width = 8, height = 6, units = "in", res = 300)
      
      par(mar = c(4, 4, 2, 1))
      
      plot(NA,
           xlim = range(hj_years),
           ylim = c(0, 1),
           xlab = "Year",
           ylab = "Observed share",
           xaxt = "n", bty = "n")
      axis(1, at = seq(2009, 2023, by = 1))
  
      for (j in seq_along(party_names)) {
        party <- party_names[j]
        pid   <- pid_list[j]
        
        y <- party_series[[party]]
        
        lines(hj_years, y, col = cols[pid], lwd = 2)
        points(hj_years, y, col = cols[pid], pch = 16)
      }
      
      legend(
        "topleft",
        legend = party_names,
        fill   = cols[pid_list],
        border = NA,
        bty    = "n"
      )
      
      dev.off()
    }
    